############################
### Placebo Test
############################

rm(list = ls())

### libraries and packages
pkg <- c("plyr", "dplyr", "tidyr", 
         "MASS", "multiwayvcov", 
         "stargazer", "lmtest", "doBy", "ggplot2", 
         "mediation", "vars", "DataCombine", "foreign",
         "mgcv", "ordinal", "reshape2", 
         "xtable", "sandwich", "rms")
lapply(pkg, require, character.only = TRUE)

# rid of E
options(scipen=999)
### setting wd for inputs
setwd("~/Dropbox/Pollution and Public Perceptions in China/Data and Analysis/Data")

# loading the data
load("ppp_cleaned")

## We conduct a placebo. We pretend that aqi.lagged is today's aqi. 

# l_sat
l_sat_p <- lm(l_sat ~ aqi.lagged + aqi.lagged2 + aqi.lagged3 + hukou + insider + soe + educ 
              + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov1 <-cluster.vcov(l_sat_p, p2$day) 
sep1 <- coeftest(l_sat_p, m.vcov1) 

# c_sat
c_sat_p <- lm(c_sat
              ~ aqi.lagged + aqi.lagged2 + aqi.lagged3 + hukou + insider + soe + educ 
              + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov25 <-cluster.vcov(c_sat_p, p2$day)  
sep2 <- coeftest(c_sat_p, m.vcov25) 

# aqi and the demand for more watchdog at the local level
l_watch_p <- lm(l_watch
                ~ aqi.lagged + aqi.lagged2 + aqi.lagged3 + hukou + insider + soe + educ 
                + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov13 <-cluster.vcov(l_watch_p, p2$day)  
sep3 <- coeftest(l_watch_p, m.vcov13) 

# aqi and the demand for more watchdog at the central level 
c_watch_p <- lm(c_watch
                ~ aqi.lagged + aqi.lagged2 + aqi.lagged3 + hukou + insider + soe + educ 
                + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov14 <-cluster.vcov(c_watch_p, p2$day)  
sep4 <- coeftest(c_watch_p, m.vcov14) # coeftest

# aqi and anti-west sentiment
anti_west_p <- lm(anti_west
                  ~ aqi.lagged + aqi.lagged2 + aqi.lagged3 + hukou + insider + soe + educ 
                  + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov25 <-cluster.vcov(anti_west_p, p2$day)  
sep6 <- coeftest(anti_west_p, m.vcov25) 

# life_sat
life_sat_p <- lm(life_sat
                 ~ aqi.lagged + aqi.lagged2 + aqi.lagged3 + hukou + insider + soe + educ 
                 + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov25 <-cluster.vcov(life_sat_p, p2$day)  
sep7 <- coeftest(life_sat_p, m.vcov25) 

# column labels
column.labels <- c("Local Government", "Central Government", "Local Oversight", "Central Oversight", 
               "Anti-West")

# covariate labels for the expanded list, 
cov.labs.c <-c("AQI", "AQI Lagged", "AQI Lagged 2","Hukou Status", "Regime Insider", "SOE Employee", "Education","Children","Married","Income", "Female", "CCP Member","Parade Period","National Holiday","Age")


# Making table with labels for full specification FOR PLACEBO
## TABLE A7
sink("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Tables/placebo.tex")
stargazer(l_sat_p, c_sat_p, l_watch_p, c_watch_p,  anti_west_p, font.size = "large", digits = 4, 
          column.labels = column.labels,
          covariate.labels = cov.labs.c,
          dep.var.labels.include = FALSE,
          model.numbers = FALSE,
          style = "qje", no.space = FALSE, keep.stat = c("n"),
          star.cutoffs = c(.10, .05, .01), 
          se = list(sep1[,2], sep2[,2], sep3[, 2], sep4[, 2],  sep6[, 2]),
          header = F, notes = "robust standard errors clustered by day in parantheses", float = F)
sink()


################################################
###### Making a figure for PLACEBO results ######
#################################################

## Figure A8
setwd("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Figures")
pdf("figure_effect2.pdf", height=8, width=6)
coef.vect <- c(sep1[,1][2], sep2[,1][2], sep3[, 1][2], sep4[, 1][2],  sep6[, 1][2])
sd.vect <- c(sep1[,2][2], sep2[,2][2], sep3[, 2][2], sep4[, 2][2],  sep6[, 2][2])
longnames <- factor(c("Local government \n satisfaction", "Central government \n satisfaction", "Demand for public oversight \n over local government",
                      "Demand for public oversight \n over central government",
                       "Belief in \n anti-west conspiracy"),
                    levels=c("Belief in \n anti-west conspiracy", 
                             "Demand for public oversight \n over central government", "Demand for public oversight \n over local government",
                             "Central government \n satisfaction", "Local government \n satisfaction"))

frame = data.frame(longnames, coef.vect, sd.vect)

c <- ggplot(frame, aes(longnames, coef.vect))
c + geom_point(size=3) + geom_errorbar(aes(ymin=coef.vect-(qnorm(.975)*sd.vect),
                                           ymax=coef.vect+(qnorm(.975)*sd.vect), width=0.2)) + coord_flip() +
  geom_hline(yintercept=0, linetype="longdash") + theme_bw() +
  xlab("Dependent Variables") + ylab("Coefficients of Daily Pollution") +
  ggtitle("Effects of Pollution on Outcomes of Interests") +
  theme(plot.title = element_text(vjust = 1.3))
dev.off()